#include <dftFFT.h>

#undef __FUNCT__
#define __FUNCT__ "DFTFFTCreate"
PetscErrorCode DFTFFTCreate(DA da, PetscReal R[], PetscReal T, DFTFFT *dft)
{
  DA             coordsDA;
  Vec            coordsVec;
  PetscReal   ***coords;
  DALocalInfo    info;
  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscTruth flag;
  dft->debug = 0;
  ierr = PetscOptionsGetInt("dft_", "-debug", &dft->debug, &flag); CHKERRQ(ierr);

  ierr = PetscObjectGetComm((PetscObject)da, &(dft->comm)); CHKERRQ(ierr);
  dft->rhoDA = da;
  ierr = DAGetInfo(dft->rhoDA, &dft->dim, 
                   dft->dims+0,dft->dims+1, dft->dims+2, 
                   PETSC_NULL, PETSC_NULL, PETSC_NULL, 
                   &(dft->numSpecies), 
                   PETSC_NULL, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr);
  if(dft->dim != 1) {
    SETERRQ1(PETSC_ERR_SUP, "Unsupported domain dimension: %d", dft->dim); 
  }
  if(dft->numSpecies <= 0 && dft->numSpecies != NUM_SPECIES) {
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number of species: %d, or NUM_SPECIES: %d", dft->numSpecies, NUM_SPECIES); CHKERRQ(ierr);
  }
  // We assume that coordinates are uniform, which is necessary for FFT
  ierr = DAGetLocalInfo(dft->rhoDA, &info); CHKERRQ(ierr);
  ierr = DAGetCoordinates(dft->rhoDA, &coordsVec); CHKERRQ(ierr);
  ierr = DAGetCoordinateDA(dft->rhoDA, &coordsDA); CHKERRQ(ierr);
  ierr = DAVecGetArray(coordsDA, coordsVec, &coords); CHKERRQ(ierr);
  // We also assume a serial DA (again, because of FFT)
  if(dft->dim == 1) {
    PetscReal *coords;
    ierr = DAVecGetArray(coordsDA, coordsVec, &coords); CHKERRQ(ierr);
    dft->L[0] = coords[dft->dims[0]-1]-coords[0];
    dft->h[0] = coords[1]-coords[0];
  }
  else if(dft->dim == 3) {
    PetscReal ***coords;
    ierr = DAVecGetArray(coordsDA, coordsVec, &coords); CHKERRQ(ierr);
    dft->L[0] = coords[dft->dims[0]-1][0][0] - coords[0][0][0];
    dft->h[0] = coords[1][0][0] - coords[0][0][0];
    dft->L[1] = coords[0][dft->dims[1]-1][0] - coords[0][0][0];
    dft->h[1] = coords[0][1][0] - coords[0][0][0];
    dft->L[2] = coords[0][0][dft->dims[2]-1] - coords[0][0][0];
    dft->h[2] = coords[0][0][1] - coords[0][0][0];
  }
  else {
    SETERRQ1(PETSC_ERR_SUP, "Unsupported domain dimension: %d", dft->dim);
  }
  ierr = DADestroy(coordsDA); CHKERRQ(ierr); 
  for(PetscInt d = 0; d < dft->dim; ++d) {
    if(dft->L[d] <= 0) {
      SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid domain size %g in dimension %d", dft->L[d], d);
    }
    if(dft->h[d] <= 0) {
      SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid mesh size %g in dimension %d", dft->h[d], d);
    }
  }
  ierr = PetscMalloc(sizeof(PetscReal)*dft->numSpecies, &dft->R); CHKERRQ(ierr);
  dft->maxR = 0;
  for(PetscInt n = 0; n < dft->numSpecies; ++n) {
    dft->R[n] = R[n];
    dft->maxR = PetscMax(dft->maxR, dft->R[n]); CHKERRQ(ierr);
  }
  if(dft->maxR <= 0) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid max species radius: %g", dft->maxR); CHKERRQ(ierr);
  }
  /* 
     Boltzmann's constant: k = 1.38e-23 [(kg m^2)/(s^2 K)] = 1.38e-5 [(kg nm^2)/(s^2 K)]
     Assuming T = 305K, we get the default kT value:
     kT = 305 [K] * 1.38e-5 [(kg nm^2)/(s^2 K)] = 4.21e-3[kg m^2/s^2] =  4.21e-21 [kg m^2/s^2]
  */
  dft->kT      = 4.21e-3;               
  dft->T = T;
  if(dft->T <= 0) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid temperature: %g", dft->T); CHKERRQ(ierr);
  }
  dft->kT = 1.38e-5*dft->T;
  

  PetscFunctionReturn(0);
}// DFTCreate()

#undef __FUNCT__
#define __FUNCT__ "DFTFFTDestroy"
PetscErrorCode DFTFFTDestroy(DFTFFT *dft) {
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = PetscFree(dft->R); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}// DFTFFTDestroy()


